/*    Copyright 2010 Tobias Marschall
 *
 *    This file is part of MoSDi.
 *
 *    MoSDi is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    MoSDi is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoSDi.  If not, see <http://www.gnu.org/licenses/>.
 */

package mosdi.util;

public class Distributions {

	private static int max3(int a, int b, int c) {
		return Math.max(a,Math.max(b, c));
	}

	private static int leftmostNonZero(double[] d) {
		int n = 0;
		while (d[n]==0.0) n++;
		return n;
	}

	/** Computes the complementary cumulative distribution, i.e. the 
	 *  right-tail sums. */
	public static double[] toCCDF(double[] distribution) {
		double[] result = new double[distribution.length];
		result[distribution.length-1] = distribution[distribution.length-1];
		for (int i=distribution.length-2; i>=0; --i) {
			result[i] = distribution[i] + result[i+1];
		}
		return result;
	}

	/** Computes the convolution of two distributions. The size of the
	 *  resulting array is chosen such that all values that are >0.0 in
	 *  double precision are contained. */
	public static double[] convolve(double[] dist1, double[] dist2) {
		// Log.printf(Log.Level.DEBUG, "convolve(%d x %d)%n", dist1.length, dist2.length);
		// Log.startTimer();
		int leftmostNonZero1 = leftmostNonZero(dist1);
		int leftmostNonZero2 = leftmostNonZero(dist2);
		double[] result = null;
		for (int n=dist1.length+dist2.length-2; n>=0; --n) {
			double d = 0.0;
			for (int i=max3(0, n-dist2.length+1,leftmostNonZero1); i<=Math.min(n-leftmostNonZero2, dist1.length-1); ++i) {
				d+=dist1[i]*dist2[n-i];
			}
			if (d>0.0) {
				if (result==null) {
					result = new double[n+1];
				}
				result[n]=d;
			}
		}
		// Log.stopTimer("convolve()");
		// TODO:
		// Log.printf(Log.Level.DEBUG, "result length: %d, leftmost non-zero: %d)%n", result.length, leftmostNonZero(result));
		return result;
	}

	/** Computes the convolution of two distributions. The size of the
	 *  resulting array equals the length of the first given distribtuion.
	 *  
	 *  @param accumulateInLastField If true, all probability mass that would leave
	 *                               the table (due to the length restriction of resulting
	 *                               array), is added to last array entry.
	 */
	public static double[] convolveLengthPreserving(double[] dist1, double[] dist2, boolean accumulateInLastField) {
		int leftmostNonZero1 = leftmostNonZero(dist1);
		int leftmostNonZero2 = leftmostNonZero(dist2);
		double[] result = new double[dist1.length];
		for (int n=0; n<result.length; ++n) {
			double d = 0.0;
			if (accumulateInLastField && (n==result.length-1)) {
				for (int i=max3(0, n-dist2.length+1,leftmostNonZero1); i<=Math.min(n-leftmostNonZero2, dist1.length-1); ++i) {
					for (int j=n-i; j<dist2.length; ++j) {
						d+=dist1[i]*dist2[j];						
					}
				}
			} else {
				for (int i=max3(0, n-dist2.length+1,leftmostNonZero1); i<=Math.min(n-leftmostNonZero2, dist1.length-1); ++i) {
					d+=dist1[i]*dist2[n-i];
				}
			}
			result[n] = d;
		}
		return result;
	}

	/** Returns the distribution given by P(X=x)=(p-1)^(x-1) * p.
	 *  Support of distribution: {1,2,...}, i.e. result[0] is always 0.0. */
	public static double[] geometricDistribution(double p, int length) {
		double[] result = new double[length];
		double d = 1;
		for (int i=1; i<length; ++i) {
			result[i] = d*p;
			d *= (1-p);
		}
		return result;
	}
}
